home *** CD-ROM | disk | FTP | other *** search
/ HPAVC / HPAVC CD-ROM.iso / pc / FAQSYS18.ZIP / FAQS.DAT / COLORQ~1.C < prev    next >
Encoding:
C/C++ Source or Header  |  1996-08-18  |  22.3 KB  |  790 lines
  1. /*
  2.  * This software is copyrighted as noted below.  It may be freely copied,
  3.  * modified, and redistributed, provided that the copyright notice is 
  4.  * preserved on all copies.
  5.  * 
  6.  * There is no warranty or other guarantee of fitness for this software,
  7.  * it is provided solely "as is".  Bug reports or fixes may be sent
  8.  * to the author, who may or may not act on them as he desires.
  9.  *
  10.  * You may not include this software in a program or other software product
  11.  * without supplying the source, or without informing the end-user that the 
  12.  * source is available for no extra charge.
  13.  *
  14.  * If you modify this software, you should include a notice giving the
  15.  * name of the person performing the modification, the date of modification,
  16.  * and the reason for such modification.
  17.  */
  18. /*
  19.  * colorquant.c
  20.  *
  21.  * Perform variance-based color quantization on a "full color" image.
  22.  * Author:    Craig Kolb
  23.  *        Department of Mathematics
  24.  *        Yale University
  25.  *        kolb@yale.edu
  26.  * Date:    Tue Aug 22 1989
  27.  * Copyright (C) 1989 Craig E. Kolb
  28.  * $Id: colorquant.c,v 3.0.1.4 1992/04/30 14:06:48 spencer Exp $
  29.  *
  30.  * $Log: colorquant.c,v $
  31.  * Revision 3.0.1.4  1992/04/30  14:06:48  spencer
  32.  * patch3: lint fixes.
  33.  *
  34.  * Revision 3.0.1.3  1991/02/06  16:46:47  spencer
  35.  * patch3: Fix round-off error in variance calculation that sometimes
  36.  * patch3: caused multiple boxes with the same center to be cut.
  37.  * patch3: Change 'fast' argument to 'flags' argument.  Currently:
  38.  * patch3:     CQ_FAST: same as old fast argument
  39.  * patch3:     CQ_QUANTIZE: data is not prequantized to bits
  40.  * patch3:     CQ_NO_RGBMAP: don't build rgbmap.
  41.  *
  42.  * Revision 3.0.1.2  90/11/29  15:18:04  spencer
  43.  * Remove a typo.
  44.  * 
  45.  * 
  46.  * Revision 3.0.1.1  90/11/19  16:59:48  spencer
  47.  * Use inv_cmap instead of find_colors.
  48.  * Changes to process multiple files into one colormap (accum_hist arg).
  49.  * Delete 'otherimages' argument -- unnecessary with faster inv_cmap code.
  50.  * 
  51.  * 
  52.  * Revision 3.0  90/08/03  15:20:11  spencer
  53.  * Establish version 3.0 base.
  54.  * 
  55.  * Revision 1.6  90/07/29  08:06:06  spencer
  56.  * If HUGE isn't defined, make it HUGE_VAL (for ansi).
  57.  * 
  58.  * Revision 1.5  90/07/26  17:25:48  rgb
  59.  * Added a parameter to colorquant for rgbmap construction.
  60.  * 
  61.  * Revision 1.4  90/07/13  14:53:31  spencer
  62.  * Get rid of ARB_ARG sh*t.
  63.  * Change a couple of vars to double so that HUGE won't cause problems.
  64.  * 
  65.  * Revision 1.3  90/06/28  21:42:56  spencer
  66.  * Declare internal functions properly.
  67.  * Delete unused global variable.
  68.  * 
  69.  * Revision 1.2  90/06/28  13:18:56  spencer
  70.  * Make internal functions static.
  71.  * Build entire RGB cube, not just those colors used by this image.  This
  72.  * is so dithering will work.
  73.  * 
  74.  * Revision 1.1  90/06/18  20:45:17  spencer
  75.  * Initial revision
  76.  * 
  77.  * Revision 1.3  89/12/03  18:27:16  craig
  78.  * Removed bogus integer casts in distance calculation in makenearest().
  79.  * 
  80.  * Revision 1.2  89/12/03  18:13:12  craig
  81.  * FindCutpoint now returns FALSE if the given box cannot be cut.  This
  82.  * to avoid overflow problems in CutBox.
  83.  * "whichbox" in GreatestVariance() is now initialized to 0.
  84.  * 
  85.  */
  86. static char rcsid[] = "$Header: /l/spencer/src/urt/lib/RCS/colorquant.c,v 3.0.1.4 1992/04/30 14:06:48 spencer Exp $";
  87.  
  88. #include <math.h>
  89. #include <stdio.h>
  90. #include "rle_config.h"
  91. #include "colorquant.h"
  92.  
  93. /* Ansi uses HUGE_VAL. */
  94. #ifndef HUGE
  95. #ifdef HUGE_VAL
  96. #define HUGE HUGE_VAL
  97. #else
  98. #define HUGE MAXFLOAT
  99. #endif
  100. #endif
  101.  
  102. static void QuantHistogram();
  103. static void BoxStats();
  104. static void UpdateFrequencies();
  105. static void ComputeRGBMap();
  106. static void SetRGBmap();
  107. #ifdef slow_color_map
  108. static void find_colors();
  109. static int  getneighbors();
  110. static int  makenearest();
  111. #else
  112. extern void inv_cmap();
  113. #endif
  114. static int  CutBoxes();
  115. static int  CutBox();
  116. static int  GreatestVariance();
  117. static int  FindCutpoint();
  118.  
  119.  
  120.  
  121. /* 
  122.  * The colorquant() routine has been tested on an Iris 4D70 workstation,
  123.  * a Sun 3/60 workstation, and (to some extent), a Macintosh.
  124.  * 
  125.  * Calls to bzero() may have to be replaced with the appropriate thing on
  126.  * your system.  (bzero(ptr, len) writes 'len' 0-bytes starting at the location
  127.  * pointed to by ptr.)
  128.  * 
  129.  * Although I've tried to avoid integer overflow problems where ever possible,
  130.  * it's likely I've missed a spot where an 'int' should really be a 'long'.
  131.  * (On the machine this was developed on, an int == long == 32 bits.)
  132.  * 
  133.  * Note that it's quite easy to optimize this code for a given value for
  134.  * 'bits'.  In addition, if you assume bits is small and
  135.  * that the total number of pixels is relatively small, there are several
  136.  * places that integer arithmetic may be substituted for floating-point.
  137.  * One such place is the loop in BoxStats -- mean and var need not necessary
  138.  * be floats.
  139.  * 
  140.  * As things stand, the maximum number of colors to which an image may
  141.  * be quantized is 256.  This limit may be overcome by changing rgbmap and
  142.  * colormap from pointers to characters to pointers to something larger.
  143.  */
  144.  
  145. /*
  146.  * Maximum number of colormap entries.  To make larger than 2^8, the rgbmap
  147.  * type will have to be changed from unsigned chars to something larger.
  148.  */
  149. #define MAXCOLORS        256
  150. /*
  151.  * Value corrersponding to full intensity in colormap.  The values placed
  152.  * in the colormap are scaled to be between zero and this number.  Note
  153.  * that anything larger than 255 is going to lead to problems, as the
  154.  * colormap is declared as an unsigned char.
  155.  */
  156. #define FULLINTENSITY        255
  157. #define MAX(x,y)    ((x) > (y) ? (x) : (y))
  158.  
  159. /*
  160.  * Readability constants.
  161.  */
  162. #define REDI        0    
  163. #define GREENI        1
  164. #define BLUEI        2    
  165. #define TRUE        1
  166. #define FALSE        0
  167.  
  168. typedef struct {
  169.     double        weightedvar;        /* weighted variance */
  170.     float        mean[3];        /* centroid */
  171.     unsigned long     weight,            /* # of pixels in box */
  172.             freq[3][MAXCOLORS];    /* Projected frequencies */
  173.     int         low[3], high[3];    /* Box extent */
  174. } Box;
  175.  
  176. static unsigned long    *Histogram,        /* image histogram */
  177.             NPixels,        /* # of pixels in an image*/
  178.             SumPixels;        /* total # of pixels */
  179. static unsigned int    Bits,            /* # significant input bits */
  180.             ColormaxI;        /* # of colors, 2^Bits */
  181. static Box        *Boxes;            /* Array of color boxes. */
  182.  
  183. /*
  184.  * Perform variance-based color quantization on a 24-bit image.
  185.  *
  186.  * Input consists of:
  187.  *    red, green, blue    Arrays of red, green and blue pixel
  188.  *                intensities stored as unsigned characters.
  189.  *                The color of the ith pixel is given
  190.  *                by red[i] green[i] and blue[i].  0 indicates
  191.  *                zero intensity, 255 full intensity.
  192.  *    pixels            The length of the red, green and blue arrays
  193.  *                in bytes, stored as an unsigned long int.
  194.  *    colormap        Points to the colormap.  The colormap
  195.  *                consists of red, green and blue arrays.
  196.  *                The red/green/blue values of the ith
  197.  *                colormap entry are given respectively by
  198.  *                colormap[0][i], colormap[1][i] and
  199.  *                colormap[2][i].  Each entry is an unsigned char.
  200.  *    colors            The number of colormap entries, stored
  201.  *                as an integer.    
  202.  *    bits            The number of significant bits in each entry
  203.  *                of the red, greena and blue arrays. An integer.
  204.  *    rgbmap            An array of unsigned chars of size (2^bits)^3.
  205.  *                This array is used to map from pixels to
  206.  *                colormap entries.  The 'prequantized' red,
  207.  *                green and blue components of a pixel
  208.  *                are used as an index into rgbmap to retrieve
  209.  *                the index which should be used into the colormap
  210.  *                to represent the pixel.  In short:
  211.  *                index = rgbmap[(((r << bits) | g) << bits) | b];
  212.  *     flags            A set of bit-flags:
  213.  *         CQ_FAST        If set, the rgbmap will be constructed
  214.  *                quickly.  If not set, the rgbmap will
  215.  *                be built more slowly, but more
  216.  *                accurately.  In most cases, fast
  217.  *                should be set, as the error introduced
  218.  *                by the approximation is usually small.
  219.  *         CQ_QUANTIZE:    If set, the data in red, green, and
  220.  *                 blue is not pre-quantized, and will be
  221.  *                 quantized "on the fly".  This slows
  222.  *                 the routine slightly.  If not set, the
  223.  *                 data has been prequantized to 'bits'
  224.  *                 significant bits.
  225.  *         
  226.  *    accum_hist        If non-zero the histogram will accumulate and 
  227.  *                reflect pixels from multiple images.
  228.  *                when 1, the histogram will be initialized and
  229.  *                summed, but not thrown away OR processed. when 
  230.  *                2 the image RGB will be added to it.  When 3 
  231.  *                Boxes are cut and a colormap and rgbmap
  232.  *                are be returned, Histogram is freed too.
  233.  *                When zero, all code is executed as per normal.
  234.  *
  235.  * colorquant returns the number of colors to which the image was
  236.  * quantized.
  237.  */
  238. #define INIT_HIST 1
  239. #define USE_HIST 2
  240. #define PROCESS_HIST 3
  241. int
  242. colorquant(red,green,blue,pixels,colormap,colors,bits,rgbmap,flags,accum_hist)
  243. unsigned char *red, *green, *blue;
  244. unsigned long pixels;
  245. unsigned char *colormap[3];
  246. int colors, bits;
  247. unsigned char *rgbmap;
  248. int flags, accum_hist;
  249. {
  250.     int    i,            /* Counter */
  251.     OutColors,            /* # of entries computed */
  252.     Colormax;            /* quantized full-intensity */ 
  253.     float    Cfactor;    /* Conversion factor */
  254.  
  255.     if (accum_hist < 0 || accum_hist > PROCESS_HIST)
  256.     fprintf(stderr, "colorquant: bad value for accum_hist\n");
  257.  
  258.     ColormaxI = 1 << bits;    /* 2 ^ Bits */
  259.     Colormax = ColormaxI - 1;
  260.     Bits = bits;
  261.     NPixels = pixels;
  262.     Cfactor = (float)FULLINTENSITY / Colormax;
  263.  
  264.     if (! accum_hist || accum_hist == INIT_HIST) {
  265.     Histogram = (unsigned long *)calloc(ColormaxI*ColormaxI*ColormaxI, 
  266.                         sizeof(long));
  267.     Boxes = (Box *)malloc(colors * sizeof(Box));
  268.     /*
  269.      * Zero-out the projected frequency arrays of the largest box.
  270.      */
  271.     bzero(Boxes->freq[0], ColormaxI * sizeof(unsigned long));
  272.     bzero(Boxes->freq[1], ColormaxI * sizeof(unsigned long));
  273.     bzero(Boxes->freq[2], ColormaxI * sizeof(unsigned long));
  274.     SumPixels = 0;
  275.     }
  276.  
  277.     SumPixels += NPixels;
  278.  
  279.     if ( accum_hist != PROCESS_HIST ) 
  280.     QuantHistogram(red, green, blue, &Boxes[0], flags&CQ_QUANTIZE);
  281.  
  282.     if ( !accum_hist || accum_hist == PROCESS_HIST) {
  283.     OutColors = CutBoxes(Boxes, colors);
  284.  
  285.     /*
  286.      * We now know the set of representative colors.  We now
  287.      * must fill in the colormap and convert the representatives
  288.      * from their 'prequantized' range to 0-FULLINTENSITY.
  289.      */
  290.     for (i = 0; i < OutColors; i++) {
  291.         colormap[0][i] =
  292.         (unsigned char)(Boxes[i].mean[REDI] * Cfactor + 0.5);
  293.         colormap[1][i] =
  294.         (unsigned char)(Boxes[i].mean[GREENI] * Cfactor + 0.5);
  295.         colormap[2][i] =
  296.         (unsigned char)(Boxes[i].mean[BLUEI] * Cfactor + 0.5);
  297.     }
  298.  
  299.     if ( !(flags & CQ_NO_RGBMAP) )
  300.         ComputeRGBMap(Boxes, OutColors, rgbmap, bits, colormap,
  301.               flags&CQ_FAST);
  302.  
  303.     free((char *)Histogram);
  304.     free((char *)Boxes);
  305.     return OutColors;    /* Return # of colormap entries */
  306.     }
  307.     return 0;
  308. }
  309.  
  310. /*
  311.  * Compute the histogram of the image as well as the projected frequency
  312.  * arrays for the first world-encompassing box.
  313.  */
  314. static void
  315. QuantHistogram(r, g, b, box, quantize)
  316. register unsigned char *r, *g, *b;
  317. Box *box;
  318. int quantize;
  319. {
  320.     register unsigned long *rf, *gf, *bf;
  321.     unsigned long i;
  322.  
  323.     rf = box->freq[0];
  324.     gf = box->freq[1];
  325.     bf = box->freq[2];
  326.  
  327.     /*
  328.      * We compute both the histogram and the proj. frequencies of
  329.      * the first box at the same time to save a pass through the
  330.      * entire image. 
  331.      */
  332.     
  333.     if ( !quantize )
  334.         for (i = 0; i < NPixels; i++) {
  335.         rf[*r]++;
  336.         gf[*g]++;
  337.         bf[*b]++;
  338.         Histogram[((((*r++)<<Bits)|(*g++))<<Bits)|(*b++)]++;
  339.         }
  340.     else
  341.     {
  342.         unsigned char rr, gg, bb;
  343.         for (i = 0; i < NPixels; i++) {
  344.         rr = (*r++)>>(8-Bits);
  345.         gg = (*g++)>>(8-Bits);
  346.         bb = (*b++)>>(8-Bits);
  347.         rf[rr]++;
  348.         gf[gg]++;
  349.         bf[bb]++;
  350.         Histogram[(((rr<<Bits)|gg)<<Bits)|bb]++;
  351.         }
  352.     }
  353.         
  354. }
  355.  
  356. /*
  357.  * Interatively cut the boxes.
  358.  */
  359. static int
  360. CutBoxes(boxes, colors) 
  361. Box    *boxes;
  362. int    colors;
  363. {
  364.     int curbox;
  365.  
  366.     boxes[0].low[REDI] = boxes[0].low[GREENI] = boxes[0].low[BLUEI] = 0;
  367.     boxes[0].high[REDI] = boxes[0].high[GREENI] =
  368.                   boxes[0].high[BLUEI] = ColormaxI;
  369.     boxes[0].weight = SumPixels;
  370.  
  371.     BoxStats(&boxes[0]);
  372.  
  373.     for (curbox = 1; curbox < colors; curbox++) {
  374.         if (CutBox(&boxes[GreatestVariance(boxes, curbox)],
  375.                &boxes[curbox]) == FALSE)
  376.                 break;
  377.     }
  378.  
  379.     return curbox;
  380. }
  381.  
  382. /*
  383.  * Return the number of the box in 'boxes' with the greatest variance.
  384.  * Restrict the search to those boxes with indices between 0 and n-1.
  385.  */
  386. static int
  387. GreatestVariance(boxes, n)
  388. Box *boxes;
  389. int n;
  390. {
  391.     register int i, whichbox = 0;
  392.     float max;
  393.  
  394.     max = -1;
  395.     for (i = 0; i < n; i++) {
  396.         if (boxes[i].weightedvar > max) {
  397.             max = boxes[i].weightedvar;
  398.             whichbox = i;
  399.         }
  400.     }
  401.     return whichbox;
  402. }
  403.  
  404. /*
  405.  * Compute mean and weighted variance of the given box.
  406.  */
  407. static void
  408. BoxStats(box)
  409. register Box *box;
  410. {
  411.     register int i, color;
  412.     unsigned long *freq;
  413.     double mean, var;
  414.  
  415.     if(box->weight == 0) {
  416.         box->weightedvar = 0;
  417.         return;
  418.     }
  419.  
  420.     box->weightedvar = 0.;
  421.     for (color = 0; color < 3; color++) {
  422.         var = mean = 0;
  423.         i = box->low[color];
  424.         freq = &box->freq[color][i];
  425.         for (; i < box->high[color]; i++, freq++) {
  426.             mean += i * *freq;
  427.             var += i*i* *freq;
  428.         }
  429.         box->mean[color] = mean / (double)box->weight;
  430.         box->weightedvar += var - box->mean[color]*box->mean[color]*
  431.                     (double)box->weight;
  432.     }
  433.     box->weightedvar /= SumPixels;
  434. }
  435.  
  436. /*
  437.  * Cut the given box.  Returns TRUE if the box could be cut, FALSE otherwise.
  438.  */
  439. static int
  440. CutBox(box, newbox)
  441. Box *box, *newbox;
  442. {
  443.     int i;
  444.     double totalvar[3];
  445.     Box newboxes[3][2];
  446.  
  447.     if (box->weightedvar == 0. || box->weight == 0)
  448.         /*
  449.          * Can't cut this box.
  450.          */
  451.         return FALSE;
  452.  
  453.     /*
  454.      * Find 'optimal' cutpoint along each of the red, green and blue
  455.      * axes.  Sum the variances of the two boxes which would result
  456.      * by making each cut and store the resultant boxes for 
  457.      * (possible) later use.
  458.      */
  459.     for (i = 0; i < 3; i++) {
  460.         if (FindCutpoint(box, i, &newboxes[i][0], &newboxes[i][1]))
  461.             totalvar[i] = newboxes[i][0].weightedvar +
  462.                 newboxes[i][1].weightedvar;
  463.         else
  464.             totalvar[i] = HUGE;
  465.     }
  466.  
  467.     /*
  468.      * Find which of the three cuts minimized the total variance
  469.      * and make that the 'real' cut.
  470.      */
  471.     if (totalvar[REDI] <= totalvar[GREENI] &&
  472.         totalvar[REDI] <= totalvar[BLUEI]) {
  473.         *box = newboxes[REDI][0];
  474.         *newbox = newboxes[REDI][1];
  475.     } else if (totalvar[GREENI] <= totalvar[REDI] &&
  476.          totalvar[GREENI] <= totalvar[BLUEI]) {
  477.         *box = newboxes[GREENI][0];
  478.         *newbox = newboxes[GREENI][1];
  479.     } else {
  480.         *box = newboxes[BLUEI][0];
  481.         *newbox = newboxes[BLUEI][1];
  482.     }
  483.  
  484.     return TRUE;
  485. }
  486.  
  487. /*
  488.  * Compute the 'optimal' cutpoint for the given box along the axis
  489.  * indcated by 'color'.  Store the boxes which result from the cut
  490.  * in newbox1 and newbox2.
  491.  */
  492. static int
  493. FindCutpoint(box, color, newbox1, newbox2)
  494. Box *box, *newbox1, *newbox2;
  495. int color;
  496. {
  497.     float u, v, max;
  498.     int i, maxindex, minindex, cutpoint;
  499.     unsigned long optweight, curweight;
  500.  
  501.     if (box->low[color] + 1 == box->high[color])
  502.         return FALSE;    /* Cannot be cut. */
  503.     minindex = (int)((box->low[color] + box->mean[color]) * 0.5);
  504.     maxindex = (int)((box->mean[color] + box->high[color]) * 0.5);
  505.  
  506.     cutpoint = minindex;
  507.     optweight = box->weight;
  508.  
  509.     curweight = 0;
  510.     for (i = box->low[color] ; i < minindex ; i++)
  511.         curweight += box->freq[color][i];
  512.     u = 0.;
  513.     max = -1;
  514.     for (i = minindex; i <= maxindex ; i++) {
  515.         curweight += box->freq[color][i];
  516.         if (curweight == box->weight)
  517.             break;
  518.         u += (float)(i * box->freq[color][i]) /
  519.                     (float)box->weight;
  520.         v = ((float)curweight / (float)(box->weight-curweight)) *
  521.                 (box->mean[color]-u)*(box->mean[color]-u);
  522.         if (v > max) {
  523.             max = v;
  524.             cutpoint = i;
  525.             optweight = curweight;
  526.         }
  527.     }
  528.     cutpoint++;
  529.     *newbox1 = *newbox2 = *box;
  530.     newbox1->weight = optweight;
  531.     newbox2->weight -= optweight;
  532.     newbox1->high[color] = cutpoint;
  533.     newbox2->low[color] = cutpoint;
  534.     UpdateFrequencies(newbox1, newbox2);
  535.     BoxStats(newbox1);
  536.     BoxStats(newbox2);
  537.  
  538.     return TRUE;    /* Found cutpoint. */
  539. }
  540.  
  541. /*
  542.  * Update projected frequency arrays for two boxes which used to be
  543.  * a single box.
  544.  */
  545. static void
  546. UpdateFrequencies(box1, box2)
  547. register Box *box1, *box2;
  548. {
  549.     register unsigned long myfreq, *h;
  550.     register int b, g, r;
  551.     int roff;
  552.  
  553.     bzero(box1->freq[0], ColormaxI * sizeof(unsigned long));
  554.     bzero(box1->freq[1], ColormaxI * sizeof(unsigned long));
  555.     bzero(box1->freq[2], ColormaxI * sizeof(unsigned long)); 
  556.  
  557.     for (r = box1->low[0]; r < box1->high[0]; r++) {
  558.         roff = r << Bits;
  559.         for (g = box1->low[1];g < box1->high[1]; g++) {
  560.             b = box1->low[2];
  561.             h = Histogram + (((roff | g) << Bits) | b);
  562.             for (; b < box1->high[2]; b++) {
  563.                 if ((myfreq = *h++) == 0)
  564.                     continue;
  565.                 box1->freq[0][r] += myfreq;
  566.                 box1->freq[1][g] += myfreq;
  567.                 box1->freq[2][b] += myfreq;
  568.                 box2->freq[0][r] -= myfreq;
  569.                 box2->freq[1][g] -= myfreq;
  570.                 box2->freq[2][b] -= myfreq;
  571.             }
  572.         }
  573.     }
  574. }
  575.  
  576. /*
  577.  * Compute RGB to colormap index map.
  578.  */
  579. static void
  580. ComputeRGBMap(boxes, colors, rgbmap, bits, colormap, fast)
  581. Box *boxes;
  582. int colors;
  583. unsigned char *rgbmap, *colormap[3];
  584. int bits, fast;
  585. {
  586.     register int i;
  587.  
  588.     if (fast) {
  589.         /*
  590.          * The centroid of each box serves as the representative
  591.          * for each color in the box.
  592.          */
  593.         for (i = 0; i < colors; i++)
  594.             SetRGBmap(i, &boxes[i], rgbmap, bits);
  595.     } else
  596.         /*
  597.          * Find the 'nearest' representative for each
  598.          * pixel.
  599.          */
  600. #ifdef slow_color_map
  601.         find_colors(boxes, colors, rgbmap, bits, colormap, 1);
  602. #else
  603.         inv_cmap(colors, colormap, bits, Histogram, rgbmap);
  604. #endif
  605. }
  606.  
  607. /*
  608.  * Make the centroid of "boxnum" serve as the representative for
  609.  * each color in the box.
  610.  */
  611. static void
  612. SetRGBmap(boxnum, box, rgbmap, bits)
  613. int boxnum;
  614. Box *box;
  615. unsigned char *rgbmap;
  616. int bits;
  617. {
  618.     register int r, g, b;
  619.     
  620.     for (r = box->low[REDI]; r < box->high[REDI]; r++) {
  621.         for (g = box->low[GREENI]; g < box->high[GREENI]; g++) {
  622.             for (b = box->low[BLUEI]; b < box->high[BLUEI]; b++) {
  623.                 rgbmap[(((r<<bits)|g)<<bits)|b]=(char)boxnum;
  624.             }
  625.         }
  626.     }
  627. }
  628.  
  629. #ifdef slow_color_map
  630. /*
  631.  * Form colormap and NearestColor arrays.
  632.  */
  633. static void
  634. find_colors(boxes, colors, rgbmap, bits, colormap, otherimages)
  635. int colors;
  636. Box *boxes;
  637. unsigned char *rgbmap;
  638. int bits;
  639. unsigned char *colormap[3];
  640. int otherimages;
  641. {
  642.     register int i;
  643.     int num, *neighbors;
  644.  
  645.     neighbors = (int *)malloc(colors * sizeof(int));
  646.  
  647.     /*
  648.      * Form map of representative (nearest) colors.
  649.      */
  650.     for (i = 0; i < colors; i++) {
  651.         /*
  652.          * Create list of candidate neighbors and
  653.          * find closest representative for each
  654.          * color in the box.
  655.          */
  656.         num = getneighbors(boxes, i, neighbors, colors, colormap);
  657.         makenearest(boxes, i, num, neighbors, rgbmap, bits, colormap, 
  658.                 otherimages);
  659.     }
  660.     free((char *)neighbors);
  661. }
  662.  
  663. /*
  664.  * In order to minimize our search for 'best representative', we form the
  665.  * 'neighbors' array.  This array contains the number of the boxes whose
  666.  * centroids *might* be used as a representative for some color in the
  667.  * current box.  We need only consider those boxes whose centroids are closer
  668.  * to one or more of the current box's corners than is the centroid of the
  669.  * current box. 'Closeness' is measured by Euclidean distance.
  670.  */
  671. static
  672. getneighbors(boxes, num, neighbors, colors, colormap)
  673. Box *boxes;
  674. int num, colors, *neighbors;
  675. unsigned char *colormap[3];
  676. {
  677.     register int i, j;
  678.     Box *bp;
  679.     float dist, LowR, LowG, LowB, HighR, HighG, HighB, ldiff, hdiff;
  680.  
  681.     bp = &boxes[num];
  682.  
  683.     ldiff = bp->low[REDI] - bp->mean[REDI];
  684.     ldiff *= ldiff;
  685.     hdiff = bp->high[REDI] - bp->mean[REDI];
  686.     hdiff *= hdiff;
  687.     dist = MAX(ldiff, hdiff);
  688.  
  689.     ldiff = bp->low[GREENI] - bp->mean[GREENI];
  690.     ldiff *= ldiff;
  691.     hdiff = bp->high[GREENI] - bp->mean[GREENI];
  692.     hdiff *= hdiff;
  693.     dist += MAX(ldiff, hdiff);
  694.  
  695.     ldiff = bp->low[BLUEI] - bp->mean[BLUEI];
  696.     ldiff *= ldiff;
  697.     hdiff = bp->high[BLUEI] - bp->mean[BLUEI];
  698.     hdiff *= hdiff;
  699.     dist += MAX(ldiff, hdiff);
  700.  
  701. #ifdef IRIS
  702.     dist = fsqrt(dist);
  703. #else
  704.     dist = (float)sqrt((double)dist);
  705. #endif
  706.  
  707.     /*
  708.      * Loop over all colors in the colormap, the ith entry of which
  709.      * corresponds to the ith box.
  710.      *
  711.      * If the centroid of a box is as close to any corner of the
  712.      * current box as is the centroid of the current box, add that
  713.      * box to the list of "neighbors" of the current box.
  714.      */
  715.     HighR = (float)bp->high[REDI] + dist;
  716.     HighG = (float)bp->high[GREENI] + dist;
  717.     HighB = (float)bp->high[BLUEI] + dist;
  718.     LowR = (float)bp->low[REDI] - dist;
  719.     LowG = (float)bp->low[GREENI] - dist;
  720.     LowB = (float)bp->low[BLUEI] - dist;
  721.     for (i = j = 0, bp = boxes; i < colors; i++, bp++) {
  722.         if (LowR <= bp->mean[REDI] && HighR >= bp->mean[REDI] &&
  723.             LowG <= bp->mean[GREENI] && HighG >= bp->mean[GREENI] &&
  724.             LowB <= bp->mean[BLUEI] && HighB >= bp->mean[BLUEI])
  725.             neighbors[j++] = i;
  726.     }
  727.  
  728.     return j;    /* Return the number of neighbors found. */
  729. }
  730.  
  731. /*
  732.  * Assign representative colors to every pixel in a given box through
  733.  * the construction of the NearestColor array.  For each color in the
  734.  * given box, we look at the list of neighbors passed to find the
  735.  * one whose centroid is closest to the current color.
  736.  */
  737. static
  738. makenearest(boxes, boxnum, nneighbors, neighbors, rgbmap, bits, colormap, 
  739.         otherimages)
  740. Box *boxes;
  741. int boxnum;
  742. int nneighbors, *neighbors, bits;
  743. unsigned char *rgbmap, *colormap[3];
  744. int otherimages;
  745. {
  746.     register int n, b, g, r;
  747.     double rdist, gdist, bdist, dist, mindist;
  748.     int which, *np;
  749.     Box *box;
  750.  
  751.     box = &boxes[boxnum];
  752.  
  753.     for (r = box->low[REDI]; r < box->high[REDI]; r++) {
  754.         for (g = box->low[GREENI]; g < box->high[GREENI]; g++) {
  755.             for (b = box->low[BLUEI]; b < box->high[BLUEI]; b++) {
  756. /*
  757.  * The following "if" is to avoid doing extra work when the RGBmap is
  758.  * not going to be used for other images.
  759.  */
  760.                     if ((!otherimages) && 
  761.                     (Histogram[(((r<<bits)|g)<<bits)|b] == 0))
  762.                     continue;
  763.  
  764.                 mindist = HUGE;
  765.                 /*
  766.                  * Find the colormap entry which is
  767.                  * closest to the current color.
  768.                  */
  769.                 np = neighbors;
  770.                 for (n = 0; n < nneighbors; n++, np++) {
  771.                     rdist = r-boxes[*np].mean[REDI];
  772.                     gdist = g-boxes[*np].mean[GREENI];
  773.                     bdist = b-boxes[*np].mean[BLUEI];
  774.                     dist = rdist*rdist + gdist*gdist + bdist*bdist;
  775.                     if (dist < mindist) {
  776.                         mindist = dist;
  777.                         which = *np; 
  778.                     }
  779.                 }
  780.                 /*
  781.                  * The colormap entry closest to this
  782.                  * color is used as a representative.
  783.                  */
  784.                 rgbmap[(((r<<bits)|g)<<bits)|b] = which;
  785.             }
  786.         }
  787.     }
  788. }
  789. #endif /* slow_color_map */
  790.